#ifndef _SPLINE_
#define _SPLINE_

#include  <iostream>
#include <cstring>     
#include  <gsl/gsl_errno.h>// header for gsl error handling
#include  <gsl/gsl_spline.h>// header for gsl splining routines 

template < class T=double>
class CSpline 
{   
 public:    
//********************************************************************
// Constructor for Spline  
CSpline(T * x_array, T * y_array, 
	int  num_pts, std::string type) :accel_ptr(NULL),spline_ptr(NULL) 
{   
  spline_type = type;    // set the private variable for the spline type
  // Allocate the accelerator  accel_ptr = gsl_interp_accel_alloc ();    // Allocate the spline according to spline_type  
  if (spline_type == "cubic ")
    {    
      spline_ptr = gsl_spline_alloc (gsl_interp_cspline, num_pts);  
    }  else
      if  (spline_type == "linear ")
	{   
	  spline_ptr = gsl_spline_alloc (gsl_interp_linear, num_pts);  
	}  
      else  {    
	std::cout << "Illegal spline type! " << std::endl;
	exit (1);   // time to quit!
      }   // Initialize the spline   
  gsl_spline_init (spline_ptr, x_array, y_array, num_pts);  
} 

~CSpline () // Destructor for Spline
{   
  // Free the accelerator and spline object
  // std::cout<<"deallocate  1 "<<std::endl;
  if(spline_ptr!=NULL )
  gsl_spline_free (spline_ptr); 
  //std::cout<<"deallocate 2  "<<std::endl; 
  if(accel_ptr!=NULL )
  gsl_interp_accel_free (accel_ptr);
} 

T y(const T x)   // find y(x)
{   
  return gsl_spline_eval (spline_ptr, x, accel_ptr);
} 

T yp (const  T x)   // find y'(x)
{   
  return gsl_spline_eval_deriv (spline_ptr, x, accel_ptr);
} 

T ypp (const  T x)   // find y''(x)
{   
  return gsl_spline_eval_deriv2 (spline_ptr, x, accel_ptr);
} 

 private:     
  std::string spline_type;       
  // type of spline 
  gsl_interp_accel *accel_ptr;    // accelerator pointer
  gsl_spline *spline_ptr;         // spline pointer
}; 

#endif
